Nonlinear parametric model for Granger causality of time series 
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We generalize a previously proposed approach for nonlinear Granger causality of time series, based 
on radial basis function. The proposed model is not constrained to be additive in variables from the 
two time series and can approximate any function of these variables, still being suitable to evaluate 
causality. Usefulness of this measure of causality is shown in a physiological example and in the 
study of the feed-back loop in a model of excitatory and inhibitory neurons. 

PACS numbers: 05.10.-a,87.10.+e,87.19.La 

I. INTRODUCTION 

Since the seminal paper by Granger [J, detecting causality relationships between two simultaneously recorded 
signals is one of the most important problems in time series analysis. Applications arise in many fields, like economy 
U, brain studies 0; S IE IE Hi > human cardiorespiratory system 0, and many others. The major approach to 
this problem examines if the prediction of one series could be improved by incorporating information of the other, as 
proposed by Granger. In particular, if the prediction error of the first time series is reduced by including measurements 
from the second time series in the regression model, then the second time series is said to have a causal influence on 
the first time series. As Granger causality was originally developed for linear systems 0, recently some attempts to 
extend this concept to the nonlinear case have been proposed. In |lOj | local linear models in reduced neighborhoods 
are considered and the average causality index, over the whole data-set, is proposed as a nonlinear measure. In 
[Tll | a radial basis function (RBF) approach has been used to model data, while in a non-parametric test of 
causality has been proposed, based on the concept of transfer entropy. A recent paper pointed out that not 
all nonlinear prediction schemes are suitable to evaluate causality between two time series, since they should be 
invariant if statistically independent variables are added to the set of input variables. This property guarantees that, 
at least asymptotically, one would be able to recognize variables without causality relationship. The purpose of 
this work is to use the theoretical results found in |13ll to find the largest class of RBF models suitable to evaluate 
causality, thus extending the results described in Moreover we show the application of causality to the analysis 
of cardiocirculatory interaction and to study the mutual influences in inhibitory and excitatory model neurons. 

Let {a;i}i=i,.,Af and {yi}i=i,.,N be two time series of N simultaneously measured quantities. In the following we 
will assume that time series are stationary. We aim at quantifying how much y is cause of x. For fc = 1 to M 
(where M = N — m, m being the order of the model), we denote a;'" — Xk+m, = {xk+m-i,Xk+m~2, ■■■,Xk), 
= {yk+m-i,yk+m-2, ■■■,yk) and wc treat these quantities as M realizations of the stochastic variables {x, X, Y) 
|14|. Let us now consider the general nonlinear model 



Wo + wi ■ * (X) + W2 • * (Y) + W3 • H (X, Y) 



(1) 



where wq is the bias term, {w} are real vectors of free parameters, $ = (tpi, (/3„^) are Ux given nonlinear real 
functions of m variables, ^ = (i/ii, are Uy other real functions of m variables, and H = (Cii ■••i Cn^„) are 
Uxy functions of 2m variables. Parameters wq and {w} must be fixed to minimize the prediction error (we assume 

e-y = TI Y.k=i (2^' - wo - wi • * (X'^) - W2 • * (Y'=) + W3 • H (X^ Y'^))' . (2) 
We also consider the model: 

x = Wo + vi • * (X) , (3) 

and the corresponding prediction error e^- If the prediction of x improves by incorporating the past values of {yi\, 
i.e. ixy is smaller than e^: then y is said to have a causal influence on x. We must require that, if Y is statistically 
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independent of x and X, then exy = at least for M — > oo. This is ensured if, for each a €! {!,..., n^y}, an index 
ia exists such that: 

Ca(X,Y) = ^,JX)r„(Y), (4) 

where Fq, is an arbitrary function of Y. As explained in PJ|, model (0, with condition is the largest class of 
nonlinear parametric models suitable to evaluate causality. We remark that e^y is equal to also at finite M , for 
statistically independent Y, if the probability distribution is replaced by the empirical measure. Exchanging the two 
time series, one may analogously study the causal influence of a; on j/. 

We choose the functions 'S' and H, in model in the frame of RBF methods thus generalizing the approach 
in We fix Ux = ny = n^y = n <C M: n centers {X'', Y''}p^]^, in the space of (X, Y) vectors, are determined by 
a clustering procedure applied to data {(X'^, Y'^)}^-^ (we use fuzzy c-means to find prototypes). We then make 
the following choice for p = 1, . . . ,n: 

^p(X) =expf-||X-XP||V2a2) , 

V'p(Y) =exp(-||Y- Y''||V2(72) , (5) 
ep(X,Y) =^p(X)^p(Y), 

a being a fixed parameter, whose order of magnitude is the average spacing between the centers. Condition Q is 
satisfied by construction. The RBF model of is recovered setting H = in (Q, i.e. it is constrained to be additive 
in variables X and Y; instead the RBF model, here proposed, can approximate any function of X and Y. 

Now we describe the application to time series of heart rate and blood pressure from patients from an intensive care 
unit, contained in the MIMIC database [igj. In healthy subjects the heart rate variability (HRV) and the systolic 
blood pressure (SBP) are interdependent. Two mechanisms determine the feed-forward influence of HRV on SBP, the 
Starling law and the diastolic decay; baroreflex regulation determines the influence of SBP on HRV 0|. We consider 
signals from six patients affected by congestive heart failure(CHF)/pulmonary edema (patients no. 212, 213, 214, 225, 
230, and 245) and six patients whose primary pathology was sepsis, a condition in which the body is fighting a severe 
infection and that can lead to shock, a reaction caused by lack of blood flow in the body (patients no. 222, 224, 269, 
291, 410, and 422). HRV and SBP time series are extracted from raw data and re-sampled at 2Hz. We fix n = 20 
and vary m from 1 to 20, taking a — 2.5/ ^/rn. Denoting x the HRV time series, and y the SBP time series, figure^] 
shows the behavior of and ey as a function of m for a typical subject. The optimal value of m corresponds to the 
knee of the curves, m — b: in terms of frequency this value corresponds to the respiratory band. In figure|21we depict 
Si ~ {ex — £xy)/£x (measuring the infiuence of SBP on HRV) and 62 — {ty — tyx)! ^y (measuring the influence of HRV 
on SBP), as a function of m, for CHF and sepsis patients. In the case of sepsis patients, the curves show a symmetric 
HRV-SBP interdependence, whilst in CHF patients the HRV SBP influence seems to be dominant, as 62 shows 
a peak at m — 5. The probability that the twelve 62 values, from all subjects and corresponding to m = 5, were 
drawn from the same population has been estimated by Wilcoxon test to be less than 10~^. The average directionality 
index D ill|, at to = 5, is equal to 0.82 for CHF patients and to —0.019 for sepsis patients. These results show, on 
one side, that sepsis condition is not characterized by unbalanced HRV-SBP loop interaction. On the other hand, 
it is well known that CHF patients show unbalanced HRV-SBP regulatory mechanism, the feed-forward HRV-SBP 
coupling being prevalent over baroreflex sensitivity ^IS] : our findings show that this effect may be described in terms 
of Granger causality between HRV and SBP time series. 

As a second application, we consider the interactions in a model of inhibitory and excitatory neurons, which has 
been studied from different points of view, and with different aims. Much attention have been dedicated to the 
mechanisms underlying modulation of visual processing by means of attention . Synchronization originating from 
feedback can be one way to explain this process of selective attention. In a previous study |2n| | , it has been shown that 
inhibitory feedback is a way for the neurons to distinguish between spatially correlated and uncorrelated input. So it 
is important to study causal influences in a closed feedback-feedforward loop, and to understand how the inhibitory 
feedback, responsible of the selectivity of the attention, modifies the causal relationships between the neurons in 
the excitatory population. Here we consider the following question: what does causality measure in the case of 
coupled firing neurons? We show that it does not merely measure coupling constants, but the combined influence 
of couplings and membrane time constants. We consider a model of interacting neurons, two excitatory and one 
inhibitory, as illustrated in flgure 13 This is a simplifled case of the models in [23 , [IJ and . All neurons are 
Leaky Integrate-and-Fire (LIF) neurons with a membrane potential V and an input current / satisfying: 

y = + (6) 

where time is measured in units of the membrane time constant and the membrane resistance is set to one. We denote 
te (tj) the membrane time constant of excitatory (inhibitory) neurons. Every time the potential of a neuron reaches 
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the threshold value Vth, a spike is fired. This resets the potential to the value Vr; it remains bound to this value for 
an absolute refractory period r^. Furthermore we normalize the reset value Vfj to zero and the threshold value Vth to 
one. The series of spikes S^{t) {j = I, 2) from excitatory neurons provide the input current // to the inhibitory LIF 
neuron, given by the convolution of the sum of the spike trains of the excitatory neurons and a standard a function, 
which mimics an excitatory post-synaptic potential (EPSP). As a consequence of the excitatory input, the inhibitory 
neuron fires action potentials. As in the case of the excitatory current, the action potential S^{t) of the inhibitory 
neuron provides the input currents Is.j to the excitatory neurons by convolution of the spike train from the inhibitory 
neuron and an inhibitory post-synaptic potential (IPSP), given by the same a function used to represent the EPSP, 
but reversed in sign. This inhibitory feedback is characterized by a gain g, and it's delivered after a delay time td. 
Summarizing, for each excitatory neuron j — 1,2 the total input current Iej is given by 

POO 

lE.jit) =fi + ilj{t) - g / dTs'{t){t -T- TD)a\e-''\ (7) 
Jo 

while the input current for the inhibitory neuron is 

Ii{t) = M + r dr (5f (i - r) + Siit ~ r)) aVe-"^ (8) 
Jo 

where is a constant base current and r\j{t) represents internal Gaussian white noise with intensity K. We use [i = 
0.5, K = 0.08, a = 18 ms and tu — 18 ms. is set to 3 ms, while te = tj — 6 ms. Note that the base current 
is below threshold, therefore the inhibitory neuron only fires in response to excitatory spike input. Moreover, using 
these parameters, the feedback model is a very stable one, resulting in firing rates from all neurons being almost 
independent on the values of membrane time constants and inhibitory feedback gain g used in this study. The firing 
rates of all the three neurons are always between 40 and 60 Hz. The model equations have been numerically solved by 
Eulero integration using a time step 0.1 ms. The causal relationships between the time series of membrane potentials 
from neurons are then evaluated using the proposed RBF approach with n = 20 and a = 2.5/ ^/m, m e [1,600]. In 
absence of inhibitory feedback, i.e. 5 = 0, there is a dominant causal influence of the excitatory input on the inhibitory 
one; on the other hand, when the feedback is turned on (g = 1.2), the causal influence is reversed, see figure 0] It is 
worth stressing, however, that in this context the measure of causality is not simply related to coupling g. Indeed, if 
we consider the same situation, but with a longer membrane time constant of the inhibitory neuron t/, the infiuence 
of the latter on the excitatory ones is notably reduced (see figure 0J. This can be explained in the following way. 
When the time constant is short, the inhibitory neuron behaves as a coincidence detector and inhibitory feedback has 
the effect of synchronizing the excitatory neurons 22], which thus tend to fire at the same time. For longer membrane 
time constant, the inhibitory neuron acts as an integrator and does not discriminate between two temporally close 
spikes; the inhibitory neuron then reacts to a smaller number of input spikes, resulting in a decreased causal infiuence 
on the two excitatory neurons. 

Another interesting situation corresponds to only one, out of the two, excitatory neuron receiving inhibitory feed- 
back: we then investigate the causal relationships between the two excitatory neurons, as mediated by the inhibitory 
one. In absence of inhibitory feedback there is no causality involved. As the feedback is switched on, as it is clear in 
figure 13 a significant causal influence, of the neuron which does not receive the feedback on the other, is observed. 
Even in this case, the causal influence depends also on the membrane time constant of the inhibitory neuron, which 
may act as a coincidence detector or as an integrator, thus considering two successive spikes from the two excitatory 
neurons as two separates spikes, or as one, respectively. This mechanism of selective feedback can be used to explain 
the phenomenon of stimulus recognition with biased attention in the visual cortex |19|. 

We have generalized the RBF approach to Granger causality so that any function of the input variables can be 
approximated. The two applications here considered show the usefulness of the proposed approach and, more generally, 
of the notion of causality. 

Acknoledgements. Valuable discussions with Nicola Ancona and Stan Gielen are warmly acknowledged. 



[1] C.W.J. Granger, Econometrica 37, 424 (1969). 
[2] J.J. Ting, Physica A 324, 285 (2003). 

[3] P. Tass et al., Phys. Rev. Lett. 81, 3291 (1998); M. Le van Quyen et al.. Brain Res. 792, 24 (1998). 

[4] S. J. Schiffet al., Phys. Rev. E 56, 6708 (1996); M. Wiesenfeldt, U. Parlitz, W. Lauterborn, Int. Jour, of Bifurcation and 

Chaos 11, 2217 (2001). 
[5] R. Sowa et al., Phys. Rev. E 71, 61926 (2005). 



4 



[6] K.J. Blinowska, R. Kus and M. Kaminski, Phys. Rev. E 70, 50902(R) (2004). 

[7] M. G. Rosenblum et al., Phys. Rev. E 64, 45202R (2001); M. Palus, A. Stefanovska, Phys. Rev. E 67, 55201R (2003). 
[8] E. Pereda, R. Quian Quiroga, J. Bhattacharya, Progress in Neurobiology 77 1 (2005). 
[9] M. G. Rosenblum et al., Phys. Rev. E 65, 41909 (2002). 
[10] Y. Chen et al., Phys. Lett. A 324, 26 (2004). 

[11] N. Ancona, D. Marinazzo and S. Stramagha, Phys. Rev. E 70, 56221 (2004). 
[12] P.P. Verdes, Phys. Rev. E 72, 26222 (2005). 

[13] N. Ancona and S. Stramaglia, cond-mat/0502511 (in press in Neural Computation). 

[14] Usually both times series are normalized in the preprocessing stage, i.e. they are linearly transformed to have zero mean 
and unit variance. 

[15] J. C. Bezdek, Pattern Recognition with Fuzzy Objective Function Algorithms (Plenum Press, New York, 1981). 
[16] http://www.physionet.org/ 

[17] Mechanisms of blood pressure waves, K. Miyakawa, C. Polosa, H.P. Koepchen (eds.). Springer, Berlin Heidelberg New York 
(1984). 

[18] CD. Pinna et al., J. Am. Coll. Cardiol. 46, 1314 (2005). 

[19] J.R. Reynolds et al., J. Neurosci. 19, 1736 (1999). 

[20] B. Doiron et al., Phys. Rev. Lett. 93(4), 048101 (2004). 

[21] D. Marinazzo, H.J. Kappen, S. Gielen, Input-driven oscillations in networks with excitatory and inhibitory neurons with 

dynamic synapses, submitted. 
[22] C. Borgers and N. KopeU, Neural Computation 15, 509 (2003). 

X 10"' 




m 



FIG. 1: ex and ty (indistinguishable) are plotted vs m, for a typical patient. 




FIG. 2; (5i (full line) and 82 (dashed line) are plotted vs m, averaged over CHF patients (above), and averaged over sepsis 
patients (below). 
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FIG. 3: Neural model architecture. 




FIG. 4: Si (full line, I ^ E) and 62 (dashed line, E I) are plotted vs m, for n — te —Gms (left), and n — 18 ms, te —6ms 
(right). Above, no feedback. Below, with feedback. 




FIG. 5; Si (full line, influence of the excitatory neuron without feedback on the one with feedback) and 82 (dashed line, influence 
of the excitatory neuron with feedback on the one without feedback) are plotted vs m, for g — 1.2. r/ — te =6ms (top); r/ = 
18 ms, TE =6ms (bottom). 



